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[1] A semianalytical algorithm for the retrieval of ice cloud properties from satellite data 
is presented. The new method is based on the semianalytical cloud retrieval algorithm and 
uses solutions of the asymptotic radiative transfer theory applicable for optically thick 
media. Therefore the new method is much less computer time expensive than the 
commonly used lookup table approaches. Usually, the cloud optical thickness and cloud 
effective droplet radius are reported for water and ice clouds even though both parameters 
are dependent on the actual crystal shape assumed in the retrieval procedures. Thus 

the authors propose to use the reduced optical thickness (ROT) and the particle absorption 
length (PAL) for the characterization of ice clouds. This implies that no a priori or 
climatological estimates of the particle shape/size distribution are necessary and increases 
the comparability of different cloud retrieval algorithms, which are built on many different 


distribution functions. If still necessary, the retrieved ROT and PAL can easily be 
transferred to values of the optical thickness and the cloud effective droplet radius by 
assuming any of those size distribution functions. The developed technique has been 
applied to data from the NASA EOS Terra Moderate Resolution Imaging 
Spectroradiometer sensor. The scene shows Hurricane Jeanne just before its landfall near 
the coast of Florida in September 2004. Both the reduced cloud optical thickness and 
the particle absorption length have been derived for the eye wall region. 


Citation: Kokhanovsky, A. A., and T. Nauss (2005), Satellite-based retrieval of ice cloud properties using a semianalytical algorithm, 


J. Geophys. Res., 110, D19206, doi:10.1029/2004JD005744. 


1. Introduction 


[2] Water exists in clouds in liquid, gaseous, and crystal- 
line forms. While remote sensing of water clouds from 
space is simplified because of the spherical shape of 
droplets and relative homogeneity of water clouds, crystal- 
line clouds are very inhomogeneous both in vertical and 
horizontal directions. In addition they consist of ice crystals 
having very complex shapes [Garrett et al., 2001; Baum et 
al., 2005a, 2005b]. This complicates retrievals of cloud 
parameters such as the size of crystals and the optical 
thickness T of a cloud. Therefore effects of inhomogeneity 
are ignored in operational ice cloud retrievals and ice 
crystals are defined in terms of some a priori assumed 
specific particle size/shape distribution functions [Platnick 
et al., 2003; Baum et al., 2005a, 2005b]. Hence a general 
approach used for liquid water cloud retrievals is closely 
followed with the substitution of the spherical particle 
model by models of crystalline media (e.g., hexagonal 
cylinders and plates). This means that the retrieved effective 
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size parameter of ice crystals depends on the assumed 
model of crystals used in the retrieval. 

[3] The task of this paper is to introduce a new param- 
eter—particle absorption length /—for ice cloud remote 
sensing applications. This parameter is habit-independent 
and can be used instead of the loosely defined ice grain size 
in the ice cloud retrieval procedures. 

[4] In order to derive ice cloud properties, the authors 
adapt the semianalytical cloud retrieval algorithm 
(SACURA) [Kokhanovsky et al., 2003] which is based on 
the asymptotic solution of the radiative transfer equations as 
described by van de Hulst [1980] to the special case of ice 
cloud remote sensing. The modified SACURA is then 
applied for the retrieval of cloud properties from Hurricane 
Jeanne using data from the Moderate Resolution Imaging 
Spectrometer (MODIS) [King et al., 1992] aboard the Terra 
and Aqua satellites. 


2. Semianalytical Cloud Retrieval Algorithm as 
Applied to Ice Clouds 
2.1. Nonabsorbing Channel 


[s] In contrast to the commonly used lookup table 
(LUT) approaches, Kokhanovsky et al. [2003] developed 
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Figure 1. Dependence of the probability of photon 
absorption 8 on the attenuation parameter c = 2aef/bo 
calculated using the Mie scattering code (circles, crosses, 
and asterisks, depending on the absorption), the Monte- 
Carlo geometrical optics code for fractals (triangles, 
depending on the absorption), and monodispersed randomly 
oriented hexagonal cylinders (hexagons, depending on the 
absorption). Lines show the dependence according to 
equation (16) at Ba = 0.47 and £ = BD.» where we find 
that B is equal to 0.9 for spheres (dashed line), 1.2 for 
hexagonal cylinders (solid line), and 1.8 for fractals (dotted 
line). The horizontal dashed line corresponds to the 
asymptotical result 8,, = 0.47. Further explanations are 
given in the text and in Table 1. 


a semianalytical cloud retrieval algorithm using the 
asymptotic solutions [van de Hulst, 1980; King, 1987; 
Kokhanovsky, 2004a] of the radiative transfer theory. 
According to van de Hulst [1980] the reflection function 
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R of an optically thick cloud in a nonabsorbing spectral 
band can be written as 


Rp, Ho, T) = Roo (uo Ho, ) — to(T) Ko(1) Kolo) (1) 


with the cosine of the solar and observation angle pọ and p, 
respectively, and the relative azimuth angle q. The diffuse 
transmittance fo(T) is given by 


1 
0.751(1—g)+a” 


(2) 


to(T) 


where g is the asymmetry parameter [van de Hulst, 1980], 
a & 1.07, T = dex L is the optical thickness, o,,, is the 
extinction coefficient, and Z is the cloud geometrical 
thickness. The escape function Ko(\1) can be described 
by the following approximate equation for h > 0.2 
[Kokhanovsky et al., 2003]: 


Kolu) = (1 + 2p). (3) 


I] we 


For a cloud over a reflective ground with Lambertian albedo 
Ago, the following equation for the cloud reflection function 
[Kokhanovsky, 2004a] applies instead of equation (1): 


R(p, Ho» o, T) = Ro (p, Hos 0) = tolT) Ko (1) Ko (Ho) 


AgotoKo(1)Ko(Mo) (4) 
1— Agoro , 


where rọ = 1 — fo is the spherical albedo for nonabsorbing 
clouds. Equations (1)—(4) reduce the calculation of the 
reflection function of a finite cloud to that of a semi-infinite 
one (R% (p, po, ©)), where R? (u, Ho, $) depends only on the 
phase function and even this dependence is rather weak 
[Kokhanovsky et al., 2003]. R&(\, o, ) itself can be 
parameterized [Kokhanovsky, 2004a, 2005] but errors 
increase with increasing satellite zenith angles and exceed 
5% for satellite zenith angles larger than 30° and solar zenith 
angles larger than 60° [Kokhanovsky, 2004a, 2005]. For nadir 
measurements the error is below 5% for all solar zenith angles 
except very low Sun positions close to the horizon which are 
not suitable for cloud retrievals anyway. While for single case 
studies and suitable observation geometries the cloud 
retrieval algorithms perform well with the approximation 
for R2.(\1, uo, 6) given by Kokhanovsky [2004a, 2005], this 
approximation is not suitable for operational full-disc 
applications using geostationary sensors. Then LUTs for the 
function R? (u, Ho, >) must be used [Nauss et al., 2005]. 


Table 1. Conditions Used in Calculations Presented in Figure 1 Obtained at n = 1.3* 


Shape Ae pm xX Method Reference 
Spheres 1-18 1.e-5, 1.e-4, 1.e-3, 1.e-2 Mie theory for polydispersed spheres Kokhanovsky [2004a] 


Fractals 8-119 
Hexagonal cylinders 45 


1.2e-5, 1.2e-4, 1.2e-3, 1.2e-2 Macke et al. [1996] 
1.2e-5, 1.e-4, 1.e-3, 5.e-2 Macke et al. [1996] 


“Spherical particles are characterized by the gamma-size distribution with the coefficient of variance equal to 0.38 [Kokhanovsky, 2004a, 2004b]. Fractals 
are presented by randomly oriented monodispersed Koch fractals of the second generation [Macke et al., 1996]. Hexagonal cylinders are randomly oriented 
and have the same size (length, 100 um; side of the hexagonal cross section, 50 um). The wavelength used in calculations was equal to 1240 nm for 
hexagonal particles and fractals. It was set equal to 550 nm for spheres. The value of aspis defined as 3v/s, where v is the average volume and s is the 
average surface of particles in the unit volume of a cloud. 


geometrical optics 
geometrical optics 
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Figure 2. Browse image of Jeanne (25 September 2004, 1615 UTC). 


Calculations of R°.(\1, uo, $) require the phase function p(0), 
which can be found using the fractal model of ice crystals as 
discussed by Macke et al. [1996] and confirmed using 
experimental measurements of crystalline media cloud phase 
functions [Kokhanovsky, 2003]. In this case, g is approxi- 
mately equal to 0.74. Very close values of g have been 
obtained by Garrett et al. [2001] from in situ airborne 
experiments. They found that g = 0.74 + 0.03 for their 
observations of crystalline media in the visible spectrum. In 
addition, Garrett et al. [2001] underlined that g is almost 
independent of the particular shape or size of ice crystals in 
this spectral region. We underline that this value of g can be 
reproduced using the fractal model of an ice crystal as 
discussed by Macke et al. [1996]. 

[6] These conclusions allow us to propose the following 
analytical equation to derive the ice cloud optical thickness 
T (see equation (2)): 


where (see equation (4)) 


No Ko(0)Ko(0o) áa J 
° R? (ðo, 9, p) — R(do, 9, p) 1 — Ag, 


or 


= ai Ko(0)K0(00) Ago | 
R? (do, Ò, p) — R(%, 9, p) 1 — Ag, 


with B, = 4/[3(1 — g)], B2 = 4 a/[3(1 — g)] leading to B, = 
5.128 and B, = 5.487 at g = 0.74. Similar analytical 
equations for the cloud optical thickness at nonabsorbing 
channels have been proposed by King [1987], who applied 
them to water cloud cases, where the value of g depends on 
the size of particles due to generally smaller radii of droplets 
as compared to dimensions of crystals. Note that it follows 
for clouds over black surfaces (Ago = 0) with account for 
equation (3): 


(1 +2u)(1 + 2p) 


+= By 7 
R? (do, 9, p) — R(90,9, p) 


Bo, (8) 


where Bf ~ 0.94. Equation (8) allows us to obtain the ice 
cloud optical thickness analytically. The results of the 
retrieval using this equation are almost insensitive to the 
size or the shape of ice crystals. 

[7] For the case of fractal particles, R2.(09, 0, p) becomes 
approximately [Kokhanovsky, 2005] 


A+ B(u + po) + Chito + p(0) 
4(H + po) 


RI (p, Hos 9) = > (9) 


where A = 1.247, B = 1.186, C = 5.157 and p(0) = 11.1 

exp(—0.0870) + 1.1 exp(—0.0140). Here 0 is the scattering 

angle in degrees found from the following equation: 0 = 

cos™'(—upo + sso cos db), where s = \/1— p, so = 
1 — pĝ. 
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Figure 3. Spectral top of atmosphere reflectance as 
calculated using SCIATRAN [Rozanov et al., 2005] for 
two modes: without gaseous absorption (thick line) and with 
gaseous absorption (dashed line). It was assumed that the 
cloud has the optical thickness 20 and is illuminated by the 
solar light at the zenith angle equal to 60°. The satellites 
observe the scene at nadir. We also accounted for the 
Rayleigh and aerosol scattering for the model atmosphere as 
described by Kokhanovsky and Rozanov [2004]. Calcula- 
tions are performed for water droplets having the gamma 
particle—size distribution with the mode radius 4 um and 
coefficient of variance 38%. 


[s] Some research groups report measured values of g 
for ice clouds as high as 0.82 [Sassen and Liou, 1979] in 
the visible range (e.g., for clouds composed of solely 
hexagonal plates). This means that the derived value of 
T is biased because of the assumption that g = 0.74. This 
may not hold for all ice clouds, especially for mixed-phase 
and multilayered clouds, so we advice to report also the 
retrieved value of the reduced optical thickness (ROT)6 = 
(1 — g) t. The reduced optical thickness is less biased 
then the optical thickness t. This is due to the fact that to 
as obtained from equation (6) is only weakly influenced by 
the shape of particles because functions R90, 0, p) and 
Ko(u) change only slightly for crystals of various forms 
[Kokhanovsky and von Hoyningen-Huene, 2004]. So the 
reduced optical thickness 6 should be mainly reported in 
outputs of the retrieval algorithms for ice cloud scenes. 
The optical thickness then can be easily estimated using 
the relationship: T = 6/(1 — g) and assumed values of g 
(e.g., g = 0.74, as discussed by Gerber et al. [2000] and 
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Figure 4. Dependence of the reflection function at A = 
1240 nm on that at A = 859 nm. 


Garrett et al. [2001]. In particular, it follows for ice clouds 
with fractal crystals: T = 46. 


2.2. Absorbing Channel 
[9] The case of an absorbing channel (e.g., at the 
wavelength A = 1.6 um) is much more complicated. 
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Figure 5. Reduced optical thickness map. See color 
version of this figure at back of this issue. 
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Figure 8. Correlation between the reduced optical thick- 
ness and cloud temperature. The cloud temperature was 
derived using the MODIS cloud retrieval algorithm. 


introduced for the nonabsorbing case. If a weakly absorbing 
wavelength is selected in a close vicinity to the visible range 
(e.g., 1.2 um), then the optical thickness, the phase function 
(and also g) are approximately equal to those parameters for a 
nonabsorbing case [Kokhanovsky, 2004a]. It means that 
(assuming that the ground albedo is known or can be 
neglected as for cloudy scenes over ocean in the near infrared 
spectral region) the value of R given by equation (15) 
depends on just one free parameter, the single scattering 
albedo wo. Therefore equation (15) can be considered as a 
single transcendent equation to find the single scattering 
albedo wo (e.g., at g = 0.74 and T given by equation (8)). 
After finding wo, one can also determine the effective 
crystal radius a, = 3v/s. Here v is the average volume of 
crystals and s is the average surface area of crystals in the 
elementary volume of a cloud. One can also introduce the 
effective diameter Def = 2a¢y: 

[10] Clearly the probability of photon absorption 8 = 1 — 
Wo increases with the size of particles. The rate of this 
increase depends on the assumption of the particular shape. 
To avoid this shape dependence, we propose to retrieve not 
the effective size of particles itself but the particle absorp- 
tion length (PAL) £ defined in such a way that 


5169) = (1 => exp(—¢/l0))B.o; (16) 
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where lo = /4nx is the so-called penetration depth 
governed by the imaginary refractive index of crystals x 
and Bo = Ploo). Note that the value of Ba for convex 
particles in random orientation coincides with those of 
spheres independently on their shape [Kokhanovsky, 2004a]. 
This allows us to derive Ba using Mie theory. In particular, 
we find Ba = 0.47 for the real part of the ice refractive index 
n = 1.3. The law given by equation (16) is characteristic for 
many types of particles in the geometrical optics scattering 
domain as found by Kokhanovsky and Macke [1997]. 
Therefore equation (15) with account for equation (16) can 
be used for the determination of the PAL Z. 

[11] The dependence off on the attenuation parameter c = 
D.y/to calculated for spherical particles, hexagonal cylin- 
ders, and fractal particles is shown in Figure 1. Details for 
the calculations are specified in Table 1. We find that all 
these diverse forms follow the common law expressed by 
equation (16) very closely and 

t= BDg, (17) 
where B depends on the habit of crystals. We find from data 
given in Figure 1: B = 0.9 for spheres, B = 1.2 for hexagonal 
prisms with parameters considered in Figure 1, and B = 1.8 
for fractal particles at n = 1.3. We found that the accuracy of 
equation (16) increases for smaller 8. Then errors of 
equation (15) are also smallest. 

[12] That means that independently on the shape assump- 
tion we can derive the spatial distribution of the particle 
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Figure 9. Spherical albedo frequency distribution. 
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Figure 11. 


4000 


shape independent retrievals of £ as proposed in this paper 


are of a certain advantage. 
3. Application of SACURA to MODIS Data 


2000 


[16] The semianalytical ice cloud retrieval algorithm 
introduced above has been applied to data from the Mod- 
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Figure 10. Frequency distribution of the top of atmo- 
absorption lengths £ for crystalline clouds with complex and 


sphere reflectance R at A = 859 nm. 


20000 


a priori unknown shapes of particles. This is a major result 


of this work. 


[13] A priori assumptions on the shape of particles can 
bias the value Dep = £/B considerably. However, this is not 
the case for £ (see Figure 1). In many cases the value of 


15000 


De itself is not needed (e.g., in climate research). Then the 
single scattering albedo is of particular importance and can 
easily be found from equation (16) using the retrieved 


value of £ (e.g., from satellite data) at a single wavelength 


(at least in the spectral region, where the real part of 
refractive index of ice n does not differ substantially from 


the value of n at the wavelength used in the retrieval of £). 


This is the case for ice in the broad spectral region (e.g., 
460-1700 nm, where n € [1.29, 1.31] [Kokhanovsky, 


2004a]). 
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[14] The value of £ is only weakly influenced by the 


imaginary part of the refractive index as it follows from 


Figure 1. However. 


the particle absorption length is gener- 


> 


ally influenced by the refraction processes on the surface of 


a scatterer. In particular. 
2/3 = 0.67 [Kokhanovs 


we have for spheres as n — 1: B 


3 


= 


0.9 at 


ky, 2004a] as compared to B 


n = 1.3. The function B(n) can be easily established for any 
particular shape or for a combination of shapes using 


geometrical optics calculations. 
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[15] However, the information on the shape distribution is 
not available, for example, in the operational ice cloud Figure 12. Statistical distribution of retrieved values of 


retrievals [Platnick et al., 2003]. Therefore we believe that 


the particle absorption length. 
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Figure 13. Dependence of the cloud optical penetration 
depth (line 1) and the cloud optical reduced penetration 
depth (line 2) on the probability of photon absorption 
[Kokhanovsky, 2004b]. 


erate Resolution Imaging Spectroradiometer (MODIS) 
[King et al., 1992] aboard of the polar orbiting Terra/Aqua 
satellites which are part of the NASA Earth Observing 
System (EOS). MODIS has 36 channels between 0.4 and 
14 um and a spatial resolution between 250 and 1000 m. 
For the present study, a MODIS scene from 25 September 
2005, 1615 UTC was used, showing Hurricane Jeanne just 
before its landfall southeast of Stuart (Florida) as a category 
3 hurricane (see Figure 2). The data were supplied by 
NASA via the Distributed Active Archive Center (DAAC, 
http://daac.gsfc.nasa.gov/) and processed within an opera- 
tional processing scheme [Nauss and Bendix, 2005]. The 
SACURA ice retrieval algorithm has been applied to the ice 
cloud regions around the eye of the hurricane which have 
been identified using NASA’s thermodynamic cloud phase 
product which is based on the brightness temperature 
difference between the 8.5 um and 11 ym MODIS channels 
[Platnick et al., 2003]. 


3.1. Reduced Optical Thickness 


[17] The 859 and 1240 nm MODIS channels have been 
selected for the cloud property retrieval over the ice cloud 
fields associated with Hurricane Jeanne (see Figure 2). It 
should be pointed out that usually the 1640 nm channel is 
used in cloud effective size retrieval procedures [Platnick et 
al., 2003]. However, ice is approximately 20 times less 
absorptive at 1240 nm then at 1640 nm, which enables the 
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sensing of larger vertical cloud columns. In addition the 
assumption wọ — | is not violated for most of the pixels. 
Therefore the simple formulation based on equation (15) 
can be used (note that both channels are almost free of 
gaseous absorption) as it can be seen from Figure 3. In 
contrast to ice clouds, the 1640 nm channel should prefer- 
ably be used for the retrieval of droplet sizes in water 
clouds. This is due to the low absorption of light by water 
droplets at 1240 nm and therefore missing information on 
the droplet size within this channel signal (see Figure 3). 
The situation is very different for crystalline media as in the 
case for the hurricane event studied in the present paper. 

[15] The dependence of the measured reflection func- 
tion at 1240 nm to that at 859 nm is presented in Figure 4 
for the central core of the hurricane. It follows that both 
reflectances almost coincide at R(859 nm) < 0.3. Then, of 
course, the retrieval of the ice crystal size using these 
channels is not possible. For thicker clouds, values of 
R(1240 nm) are lower than those at 859 nm because of 
the increased importance of light absorption effects for 
larger wavelengths by large ice crystals. In particular, light 
absorption at 859 nm by ice is two orders of magnitude 
lower than that at 1240 nm. So the level of deviations of 
points from a straight line in Figure 4 can be used to 
estimate the size of crystals as discussed in the previous 
section. The optical thickness is determined from mea- 
surements at 859 nm using equation (5). We suppose that 
the value of T derived at 859 nm coincides with that at 
1240 nm. This is an accurate assumption due to the 
spectral neutrality of the ice cloud extinction coefficient 
in the near IR. 

[19] A statistically insignificant number of points deviate 
from the general rule R(859) > R(1240 nm) (see Figure 4). 
These points are not consistent with the model of a thick 
plane-parallel cloud over a black surface under the assump- 
tion of very low levels of gaseous absorption used in this 
paper. So we conclude that either our model is not valid for 
these selected pixels of a hurricane or there are problems 
with the measurements. Most likely, the deviation is caused 
by 3-D radiation effects especially in the eye wall region 
that are not accounted for in our retrieval model. The 
derived map of the ROT is given in Figure 5. We see that 
values of 6 are larger than 15 in the vicinity of the core. 
They are smaller than 10 on the periphery. The correspon- 
dent frequency histogram is shown in Figure 6. It follows 
that the frequency distribution /(6) has two maxima, reflect- 
ing different conditions in the central parts of the hurricane 
with the modal value of 6 = 14 as compared to the 
periphery where the modal value of 6 is close to 3. 
Retrievals for 6 < 1 are less accurate because of the poor 
accuracy of approximation (4) for thin clouds. 

[20] The correlation plot for the retrieved reduced 
optical thickness against the MOD06 algorithm derived 
cloud optical thickness T,, is presented in Figure 7. One 
must expect that 6 = (1 — g)t,,. This is indeed the case 
as one finds 7,, = 50 from data given in Figure 7, which 
corresponds to g = 0.8. 

[21] We do not show the SACURA-derived spatial dis- 
tribution of the optical thickness t because it is influenced 
by the vertical distribution of the parameter e = 1 — g, 
which presumably changes from 0.25 for ice clouds in the 
upper part of the hurricane to 0.15 for water clouds 
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Figure 14. Statistical distribution of retrieved values of (a) the single scattering albedo and (b) the 
frequency distribution of the top of atmosphere reflectance R at A = 1240 nm. 


underneath. Depending on assumptions of the vertical 
profile g(z), where z is the vertical coordinate, considerably 
different values of T can be derived. 

[22] The correlation of 6 with the MODIS-derived cloud 
temperature T is shown in Figure 8. Even there is no such 
clear correlation as in Figure 7, a tendency exists for thicker 
clouds to have smaller cloud top temperatures. This can be 
easily understood because of the fact that such clouds are 
also geometrically thicker and therefore they penetrate to 
higher altitudes, where Tis rather small. The minimal values 
of Tare close to —80°C in the case studied. It follows that T< 
—20°C for most of retrievals confirming that we really 
consider ice and not water clouds in this study. 

[23] Equation (6) allows us to find the spherical albedo of 
a hurricane rọ = 1 — tọ at the wavelength 859 nm. This is 
shown in Figure 9. It follows that the modal value of ro is 
close to 0.91 for the hurricane core. The frequency distri- 
bution Aro) increases rather slowly with ro and then abruptly 
drops out at rọ œ% 0.92. This means that the frequency 
distribution presented in Figure 9 follows a beta distribution 
law common for bounded processes (e.g., the maximal 
value of ro is 1 by definition). Similar trend is seen in the 
reflection function distribution presented in Figure 10. 


3.2. Particle Absorption Length 

[24] We do not present the effective crystal size in this 
section, which is acommon practice for modern cloud remote 
sensing techniques [Platnick et al., 2003], but rather the 
particle absorption length, which is much more robust against 
assumptions of crystals habits. To find the PAL, we need to 


know not only 6 but also T (see equations (11) and (14)). 
The transfer from 6 to t makes no problem for pure 
crystalline clouds (e = 0.25) but in the case considered 
here we face a multilayered cloud system with both ice 
crystals and water droplets. So the retrieval generally 
becomes much more complicated and, as it was specified 
above, the profile e(z) is not known a priori. To avoid 
this problem, we assumed in equation (11): € = 1 — g = 
0.25. This may bias results for thin clouds but actually 
has no effect on most of our retrievals because for the 
hurricane case the cloud thickness is quite large and the 
asymptotic limit of a semi-infinite cloud is almost reached 
at the 1240 nm. Therefore the value of T is not a crucial 
parameter of the retrieval in this specific case. This also 
means that the two last terms in equation (11) make only 
a small contribution for optically thick clouds with large 
crystals. Interestingly, equation (15) can be solved ana- 
lytically with respect to y (and therefore 8, £) as T — oo. 
Then the last term in equation (15) can be dropped. 

[25] The map of the particle absorption length derived 
from our bispectral algorithm as underlined above is 
shown in Figure 11. We see that in most of the cases 
the PAL is in the region 100-150 jm. The spatial 
distribution of £ does not have the symmetry with respect 
to the core. In particular, larger values of £ and maybe 
also larger hurricane momentum and larger mass, are 
observed in the western part, where the direction of the 
hurricane propagation is located. It could be of a practical 
value to confirm larger values of the PAL in the hurricane 
propagation direction when studying more hurricane 
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length and the effective radius a, The values of aep are 
derived using the MODIS cloud retrieval algorithm. 
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events. Note that clouds inside the eye are thin and 
therefore retrievals of both / and 6 using SACURA are 
not reliable. Then the standard lookup table approach 
must be applied. 

[26] The correspondent frequency histogram is shown in 
Figure 12. It follows that the frequency distribution A£) has 
the maximum at £ ~ 0.14 mm and the distribution A£) has 
just one mode as compared to f(6) (see Figure 6). This 
confirms a more uniform distribution of crystal sizes in the 
hurricane as compared to the total water distribution. It 
should be remembered that the wavelength 1240 nm used 
for the particle sizing has a finite penetration depth. So our 
estimations are limited to depths L < L,, where L, is the 
geometrical penetration depth defined as Lp = T,/0 ex, Where 
O ex, 18 the cloud extinction coefficient and 7, is the optical 
penetration depth defined by Kokhanovsky [2004b] as the 
value of optical thickness such that R(T) = 0.9R(00), where 
R(oo) is the reflection function for a semi-infinite layer. 
Clearly, t, depends on wọ. In particular, we present the 
dependence of T, (and also the reduced penetration optical 
thickness 6, = (1 — g)t, in Figure 13 for nadir observation 
and the solar zenith angle of 60°. It was assumed that g = 
0.75 and the phase function has been approximated by that 
for fractal particles [Macke et al., 1996]. The frequency 
distribution of wo derived using SACURA is shown in 
Figure 14a. The distribution of the single scattering albedo 
has a similar shape as the frequency distribution of the 
reflection function R(1240 nm) shown in Figure 14b. 

[27] It follows from Figure 14a that the modal value of 
the single scattering albedo is equal to 0.9925. This con- 
firms the applicability of the weak absorption approxima- 
tion to the case studied and also allows the estimation of 7, 
using equation (8) given by Kokhanovsky [2004b]. One 
obtains from Figures 13 and 14a that the modal value of 6 
equals 0.0075, the modal value of T, is approximately equal 


0.20mm 


Figure 16. 
2004). 


Images of ice crystals in hurricane Humberto. 2001 (A. Heymsfield, private communication, 
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Figure 17. Correlation between the particle absorption 
length and the reduced optical thickness as derived by 
SACURA. 


to 10, and the modal value of ô, approximately equals 2.5 at 
» = 1240 nm. 

[28] Taking into account that the modal value of 6 at \ = 
859 nm equals 14 for the hurricane core, we obtain that the 
wavelength 1240 nm senses approximately 20% of the 
hurricane vertical extent starting from its top (6, ~ 2.5). It 
means that remote sensing of crystals at the wavelength 
1240 nm allows us to perform an averaging for quite large 
vertical columns. 

[29] Because of the general increase in the size of the ice 
crystals with the distance from the cloud top height down- 
ward, one may expect that sizes derived for less penetrated 
wavelengths are generally smaller. In particular, we show 
the correlation plot between the effective crystal radius as 
derived by the MOD06 product at 1640 nm versus the 
SACURA-derived PAL at A = 1240 nm in Figure 15. We 
see a clear correlation between these two parameters with 
values of £  10a.y: 

[30] Unfortunately, we are not aware of any aircraft 
campaigns to collect particles in the hurricane studied here. 
However, such results are available for Hurricane Humberto 
(see Figure 16). The average maximum dimension of 
crystals collected at the top of Humberto (0.15 mm) is 
close to that given in Figure 12. So we may assume that the 
PAL may well correlate with the maximal crystal dimen- 
sions Dmax measured in the field campaigns. This assump- 
tion must be checked in future research. 


KOKHANOVSKY AND NAUSS: CLOUD RETRIEVAL ALGORITHM 


D19206 


[31] In conclusion, we present correlation plots of the 
derived values of £ with 6 and T (see Figures 17 and 18). It 
follows that the particle absorption length decreases with the 
reduced optical thickness and the value of £ increases with 
temperature. Both facts indicate that crystals at higher 
atmosphere levels have smaller sizes. In particular, it 
follows from Figure 18 that 


l=ae", (18) 


where a = 550 um, b = 0.02 and T is given in °C. A similar 
expression (but for Dmax) has been found by Ivanova et al. 
[2004] using a limited number of in situ measurements. 
They reported the following values of parameters: a = 
332.58 um and b = 0.0335. An interesting consequence of 
equation (18) is the possibility to estimate (in the statistical 
sense) the cloud temperature from the measurements of the 
particle absorption length. 


4. Conclusions 

[32] The semianalytical algorithm proposed by 
Kokhanovsky et al. [2003] has been extended to the case of 
ice clouds. Our retrieval technique is applicable only to the 


case of optically thick extended cloudiness as selected for this 
study on the basis of MODIS measurements. The paper gives 
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Figure 18. Correlation between the particle absorption 


length and the cloud temperature. The cloud temperature 
was derived using the MODIS cloud retrieval algorithm. 


11 of 13 


D19206 


b) 0.6 


0.5 


0.4 


reflection function 


E 
2 
= 

© 

< 

E 
2 

E 
2 
= 

O 

O 
= 

2 


0.1 


optical thickness 


Figure Al. 


KOKHANOVSKY AND NAUSS: CLOUD RETRIEVAL ALGORITHM 


10 


optical thickness 


D19206 


c) os 


reflection function 


0.1 


0.0 


100 10 


optical thickness 


100 


Dependence of reflection function on the optical thickness calculated using the exact 


radiative transfer code for the black underlying surface, nadir observation, different solar angles (30°, 
45°, 60°, and 75°) and (a) wo = 1, (b) wo = 0.99, and (c) wo = 0.98 for the phase function of fractal ice 
crystals as described by Mishchenko et al. [1999]. Dotted lines give results as obtained from the 
analytical equation (15) with account for equations (12)—(14). The value of g is equal to 0.7524 
[Mishchenko et al., 1999] for the phase function under study. 


the first application of SACURA to ice cloud fields. The 
results obtained allow us to study the complex structure of the 
hurricane with a spatial resolution of 1 km not accessible to 
microwave techniques usually used to monitor hurricanes 
(see Figures 5 and 11). 

[33] We proposed two new parameters to be remotely 
sensed from satellites over thick ice clouds: the particle 
absorption length and the reduced optical thickness. These 
parameters are less biased as compared to aep and T and 
therefore allow us to avoid the great uncertainty associated 
with the retrieval of crystal sizes from airborne and space- 
borne optical instruments. We underline that the retrieved 
values of a¿y and T strongly depend on a priori assumptions 
with respect to the habits of crystals. 

[34] First maps of the spatial distribution of the reduced 
optical thickness and the particle absorption lengths are 
derived using MODIS bispectral measurements at 859 nm 
and 1240 nm free of gaseous absorption. For the case studied, 
the reduced optical thickness has two modes (6 ~ 3, 6 14) 
and the most frequent value of the particle absorption length is 
0.14 mm. 

[35] The information on typical sizes of ice crystals in 
hurricanes is of importance in understanding the hurricane 
dynamics and also for weather prediction, radiative transfer 
models, convective dehydration, and stratospheric chemis- 
try [Webster and Heymsfield, 2003]. Large values of ¢ 
obtained in our work suggest large near-surface condensate 
transported to upper levels by deep convection. Normally, 
crystal sizes are much smaller for ice clouds originated close 
to the tropopause [Garrett et al., 2004] because of many 
factors including first of all the shortage of water vapor at 
high altitudes. 

[36] The retrieved values of £ for thinner clouds may be 
biased because of the possible broken clouds conditions. At 
this stage it is difficult to relate the value of @ to the 


error, % 


10 


15 20 25 


optical thickness 


Figure A2. Errors of equation (15) calculated under the 
same conditions as in Figure Al. The first number in 
parentheses gives the observation angle in degrees, and the 
second number is the single scattering albedo. 
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geometrical characteristics of particles. However, we may 
assume that £ = Dmax, Where Dmax is the maximum 
dimension of the particle. This is approximately the case 
for ice spheres, where £ ~ 0.9Dmax (see Figure 1). Here 
Dmax is equal to the average diameter for spheres. 

[37] The retrieved values of / must be in the range of 
experimentally measured values of Dmax derived from in 
situ experiments. We found that the most frequently occur- 
ring values of £ derived using SACURA well correspond to 
in situ measurements of the maximal size of crystals in a 
hurricane (see Figure 16 and also http://earthobservatory. 
nasa.gov/Newsroom/NasaNews/2002/200204299301.html). 
Clearly, this assumption should be validated in future. The 
modal values of Dmax are usually in the range 50-300 pm 
as was found by in situ measurements [Baum et al., 2005a]. 
The majority of the retrieved particle absorption lengths £ 
also belongs to this interval (see Figure 12). 


Appendix A: Accuracy of Approximate Equation 
for the Reflection Function 


[38] The accuracy of main equation (15) is studied in 
Figures Al and A2 for a specific case of the black 
underlying surface, nadir observation and wavelengths 
859 and 1240 nm for the several values of wo and 
selected solar zenith angles. The fractal model of scatter- 
ers as discussed by Macke et al. [1996] and Mishchenko 
et al. [1999] has been used to calculate the phase 
function. We find that the error increases for smaller 
values of the single scattering albedo and the optical 
thickness, being smaller than 5% at T > 6 and wo > 
0.98. The values of T and also wo at 1240 nm were larger 
than values reported above for a great majority of cases 
studied in this work. Therefore we estimate that the error 
of equation (15) was below 2% for most of retrievals 
performed. Even smaller errors (typically, below 1%) 
occur for the parameter 8 calculated using equation (16) 
as compared to the ray tracing code at wo > 0.98. The 
parameter g does not change considerably in the spectral 
range from 859 to 1240 nm for fractal particles and close 
to the assumed value of 0.74 (within 1%). Therefore we 
conclude that our errors are generally below or close to 
calibration errors of MODIS. 

[39] The error propagation analysis for SACURA has 
been performed by Kokhanovsky et al. [2003]. Therefore 
we do not repeat this analysis here and refer to our previous 
work. 
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Figure 5. Reduced optical thickness map. Figure 11. Retrieved particle absorption length map. 
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